################################################################################################

options(stringsAsFactors=F)
library(optparse)
library(parallel)
library(data.table)
library(Seurat)
library(Rmagic)

## 很关键
Sys.setenv(RETICULATE_PYTHON = "~/miniconda3/envs/magic_v1/bin/python3")
library(reticulate)
# repl_python()

##########################################################################################

option_list <- list(
    make_option(c("--cluster"), type = "character"),
    make_option(c("--magic_dir"), type = "character"),
    make_option(c("--output_dir"), type = "character")
)

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

clustern <- gsub( "cluster" , "" , opt$cluster)
magic_dir <- opt$magic_dir
output_dir <- opt$output_dir

################################################################################################

load(paste0(magic_dir,"/cluster",clustern,"_magic_for_correlation.Rdata"))

################################################################################################

cluster_cor <- data.frame(rbindlist(mclapply(1:nrow(cluster_overlap_effc),function(nrowi){
	# print(nrowi)
	atac_id <- cluster_overlap_effc$V4[nrowi]
	rna_id <- cluster_overlap_effc$V11[nrowi]
	atac_exp <- as.numeric(atac_data[atac_id,])
	rna_exp <- as.numeric(rna_data[rna_id,])
	cor_res <- cor.test(rna_exp,atac_exp,method = 'spearman')
	return(data.frame(cluster_overlap_effc[nrowi,],correlation=cor_res$estimate,pvalue=cor_res$p.value))
},mc.cores=10)))
cluster_cor$FDR <- p.adjust(cluster_cor$pvalue, method = "fdr")
cluster_cor$classt <- ifelse(cluster_cor$correlation>0&cluster_cor$correlation<0.1,"0-0.1",ifelse(cluster_cor$correlation>0.1&cluster_cor$correlation<0.2,"0.1-0.2",ifelse(cluster_cor$correlation>0.2&cluster_cor$correlation<0.3,"0.2-0.3",ifelse(cluster_cor$correlation>0.3&cluster_cor$correlation<0.4,"0.3-0.4",ifelse(cluster_cor$correlation>0.4&cluster_cor$correlation<0.5,"0.4-0.5",ifelse(cluster_cor$correlation>0.5&cluster_cor$correlation<0.6,"0.5-0.6",ifelse(cluster_cor$correlation>0.6,">0.6","<0")))))))
# table(spg_cor1$V13,spg_cor1$classt)
# spg_cor1$cor_type <- ifelse(spg_cor1$pvalue<0.05&spg_cor1$correlation>0,"positive",ifelse(spg_cor1$pvalue<0.05&spg_cor1$correlation<0,"negative","none"))
save(cluster_cor,file=paste0(output_dir,"/cluster",clustern,"_correlation_magic.Rdata"))

